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NOMENCLATURE 

A Finite difference coefficient 

C Convective flux 

D Diffusive flux 

d Cavity depth 

F Right hand side of the equation L<j> = F 

I Restriction/prolongation operator 

p Pressure 

R Residual in the solution 

Re Reynolds number (u w d/v) 

S Source term 

u x-direction velocity 

u w Top wall u velocity 

v y-direction velocity 

x,y Coordinate directions 

p Fluid density 

e Error tolerance 

n Smoothing factor 

v Fluid kinematic viscosity 

<j> General transport variable 
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STUDY OF SECOND ORDER UPWIND DIFFERENCING 
IN A RECIRCULATING FLOW 


by 

S. P. Vanka 
ABSTRACT 

In this study, we have investigated the accuracy and stability 
of the second order upwind differencing scheme originally suggested 
by Price et al. , and recently used by several others. The solution 
algorithm employed is based on a coupled solution of the nonlinear 
finite-difference equations by the multigrid technique. 
Calculations have been made of the driven cavity flow for several 
Reynolds numbers and finite-difference grids. We observe that in 
comparison with the hybrid differencing, the second order upwind 
differencing is somewhat more accurate but it is not monotonically 
accurate with mesh refinement. Also, the convergence of the 
solution algorithm deteriorates with the use of the second order 
upwind differencing. 


1. INTRODUCTION 

The calculation of practical fluid flows through the numerical solution 
of the governing partial differential equations relies heavily on the 
accuracy and stability of the finite difference (or finite element) 
scheme. In several earlier studies [1-22], existing schemes have been 
assessed and modifications have been suggested. However, the quest for an 
optimal discretization procedure still continues. The main difficulty 
arises from the (nonlinear) convective terms expressed by first order 
spatial derivatives of the flow variables. For low cell Reynolds numbers 
(less than two), the central difference operator is both stable and second- 
order accurate. However, in many practical calculations, the cell Reynolds 
numbers are much larger than this value. If central differencing is used 
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when cell Reynolds numbers are greater than two, spurious oscillations are 
known to result. The suggested cure [17] is to use a blend of central and 
upwind differences, giving the acronymed "hybrid differencing scheme.” The 
hybrid differencing scheme has been used in several earlier studies but has 
been shown to contain significant numerical diffusion. Unless considerable 
mesh refinement is undertaken, the solution can be in appreciable error. 

The second order upwind differencing (original idea traced to Price et 
al. , [23]) has been suggested as an alternate to first order upwind 
differencing. In this scheme (for the convective term) the value at a local 
node is connected to two upstream values rather than just one value. This 
scheme can be formally shown to have second order accuracy. Atias et al., 
[8] have used this scheme to study the driven cavity flow and the problem of 
impinging jet on a normal flat plate in conjunction with the stream-function 
vorticity approach. Wilkes and Thompson [24] have used the second order 
upwind scheme in the calculation of laminar and turbulent flow in sudden 
expansions and contractions. Recently Shyy and Correa [1] have used the 
scheme for the laminar driven cavity problem and other more complicated 
situations. These three studies have concluded that the second order upwind 
scheme is more accurate than the first order upwind formulation and can 
therefore be a substitute for the hybrid scheme. 

In two parallel developments, the Quadratic Upstream Differencing 
Scheme (QUICK) and the Skew Upwind Differencing Scheme (SUDS) have been 
suggested as alternates to the first order differencing. In QUICK [21], the 
concept is to evaluate the convective flux out of an interface of a finite 
difference cell by fitting a parabola between two upstream and one 
downstream nodes containing that interface. In SUDS [22], on the other 
hand, the concept is different. Here first order upwind differencing is 
used, but it is applied along the skewed streamline passing through the 
interface in question. Thus, upwinding is used in a vector sense, rather 
than along the resolved flow directions. Both QUICK and SUDS have been 
evaluated In a number of recent studies [1,5,9] and are shown to be more 
accurate than the first order upwind scheme. However, in using these 
schemes, difficulties of convergence and overshooting of the results from 
known exact solutions have been reported. Recently [5] modifications to 
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SUDS scheme in the form of bounding the solutions have been investigated 
with moderate success. 

In this study, the second order upwind differencing is incorporated 
into a recently developed efficient solution algorithm for computation of 
general fluid flows. The new solution algorithm [25,26] is significantly 
different from many others such as the SIMPLE (or its variants) algorithm 
used in several earlier investigations, in the form of the TEACH 
computer program [27]. The present algorithm, BLIMM (for Block- Implicit 
Multigrid Method), is based on a coupled solution of the momentum and 
continuity equations by the powerful multigrid technique. Used in 
conjunction with the hybrid differencing scheme (HDS) it was observed to 
provide much faster convergence than SIMPLE, converging typically in ten to 
twenty iterations even with finite difference grids as large as 512x512 
nodes in the square cavity problem. Consequently, the CPU times were 
observed to be much smaller than those of the TEACH computer program. In 
this study, the stability and accuracy of the higher order scheme are 
critically assessed by comparing its performance with the hybrid scheme. 
The calculations have been made for the flow in a two dimensional square 
cavity with a moving top wall. 

The present paper is outlined as follows. First, the second order 
upwind differencing is detailed and contrasted with the hybrid finite 
differencing scheme. Three alternate formulations of second order upwind 
differencing are shown. In Section 3, the solution algorithm is briefly 
explained. Section 4 details the calculations made for the driven cavity 
problem and their results. The experiences are discussed in Section 5. 


2. SECOND ORDER UPWIND DIFFERENCING 


Consider the one-dimensional scalar transport equation 


u4> 

X 


v* = 0 

XX 


9 


( 1 ) 
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with some given boundary conditions. The finite difference form of the 
diffusion term can be expressed without any difficulty by the central 
difference operator. Thus 


"fr u i+i ' 2 *i + W • 


5x 


(2) 


where i is the index counter for the finite difference nodes in x 
direction. 6x is the mesh size. 


Similarly, the convective term (u$ x ) can be expressed by a central 
difference operator 


^i = 


25x (<|, i+l 


- W 


(3) 


However, with this formulation when the cell Peclet number (u^x/v) is 
greater than 2, the coefficient of <j>^ becomes negative. Consequently, the 
upwind formulation 


= &r ( *i " W 5 u i > 0 » 


-■& ( * i+ i " V ■ u i < 0 


(4) 


is frequently used. The above upwind formulation has only first order 
accuracy compared with the second order accuracy for central differencing. 


The second order upwind formulation which has a higher order truncation 
error is as follows 

<“♦*>! (3 *1 ' 4 +i-l + W ’ “i > ° 

' 2& <_3<J, ± + 4 *1+1 ‘ *i+2 > ’ U i 4 ° ' <5> 

This scheme has formally 0(5x 2 ) truncation error. Wilkes and Thompson [24] 
have shown that the scheme has quadratic convergence to the exact solution, 
compared with linear convergence of the first order upwind scheme. 
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In order to implement the second order upwind scheme, a deferred 
correction procedure is frequently followed. That is, part of the 
differencing is written in terms of the coefficients linking the neighbor 
values and the remaining is treated explicitly as a source term. This 

enables a formulation that has a triagonal (in one dimension) or a 
pentadiagonal (two dimensions) matrix structure which is easy to solve. The 
deferred correction procedure is similar to the "curvature correction” used 
by Han et al. , [13] in conjunction with the QUICK scheme. 


The second order upwind differencing can be implemented in at least 
three ways, with varying degrees of implicitness and corresponding source 
terms. They are: 


Scheme 1 


= &T ( *i ‘ W + u i 


— i -^-+ 0( Sx; for u.> 0 , 

26x l 


’ £ ( *i + l ' *i> - -k ( *l + 2 ' 2 *i + l + *i> + 0(5x2) £ ° r V 0 ' 

( 6 ) 


The first term is the upwind formulation; the second term is the explicit 
correction. The advantage with this formulation is that it can be easily 
implemented in a computer program that already uses the upwind differencing. 

Scheme 2 


Scheme 2 has been used by Wilkes and Thompson [24] and the terms are 
arranged as follows. 


^x^i 2dx ( *i *i-l^ + 2Sx ( *i-2 

= 2& ( *i+l *i ) 2&T ( ^i+2 


~ 4> 4 ,) + 0(5x 2 ) for u.> 0 , 
i-1 i 

- 4> i+ ^ ) + 0(5x 2 ) for u i < 0 . 


(7) 


In this formulation more implicitness is brought into the coefficient terms. 
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Scheme 3 


Scheme 3 has been used by Atias et al. , [8] and Shyy and Correa [1] 

with the following rearrangement. 





fix 

2u. 

fix 


^i " ^ + 2 fix (<1) i-2 

u i 

( *i+l " *i* 2 fix (<t> i+2 


- $ ) + O(fix^) for u ± > 0 , 

2 

- <^) + 0(5x ) for u^< 0 . 


( 8 ) 


Scheme 3 has the largest coefficient term, having values twice those for 
Scheme 1. The three different formulations offer different degrees of 
implicitness in the solution of the equations. Thus their stability 
characteristics can be different. However, the final solution to the 
finite-difference equations is the same. 


A. Boundary Conditions 

Because the second order upwind differencing connects two upstream 
nodes, it is necessary to follow special practices for the finite difference 
nodes adjacent to a boundary; otherwise reference will have to be made to a 
value outside the flow domain. Atias et al., [8] have recommended switching 
back to first order upwind differencing for cells adjacent to a solid 
boundary. In Scheme 1 this implies setting the source term to zero for 
finite difference nodes i = 2 and i = (IMAX-1) (for a one dimensional 
problem). Shyy and Correa [1] have used the hybrid formulation near the 
boundaries. The hybrid formulation at boundaries is the same as the upwind 
formulation except that when the boundary cell Peclet number is less than 
two, central differencing is used. In the present work, the upwind 
differencing scheme, as used by Atias et al. , is employed for the near- 
boundary cells. 


B. Present Implementation 

In this study, the second order upwind differencing is evaluated by 
calculating the laminar flow field in a square cavity with a moving top 
wall. The equations governing the flow are the two-dimensional momentum and 
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mass continuity equation for a planar flow, expressed as follows. 


(uu) + (vu) = + v(u + u ) , 

x y p8x xx yy 

(vu) + (w) = + v(v + v ) , 

x y p8y xx yy' 


u + v = 0 . 
x y 


( 9 ) 

( 10 ) 

(ID 


A staggered mesh system has been used for storing the velocities and the 
pressures. The finite-difference equations are derived by integrating them 
over control volumes surrounding the velocities and the pressure. Because 
of the staggered mesh arrangement, these control volumes are different for 
the three equations, as shown in Figure 1. 


In the present study, we have considered only Schemes 1 and 3. The 
performance of Scheme 2 is expected to be in between that of Scheme 1 and 
Scheme 3. In Scheme 1, the coefficients are assembled exactly as those for 
upwind differencing. The source term is calculated from existing values in 
store and is added to the right hand side. The general form of the finite 
difference equation is 


C ^i , j “N*i,j+1 


where <j>. . represents u. . and v. .. For a uniform mesh of 6x and 5y 
intervals in x and y directions, A^j, A g , etc. are given by 


A c = a n + A s + \ + \ ’ 

(13) 

A w = max(0, C x _) + D x _ , 

(14) 

A e = max(0, - C x+ ) + D x+ , 

(15) 

A = max( 0, C ) + D , 

(16) 

S y- y- 


Ajq - raax (0, - C y+ ) + D y+ , 

(17) 
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where 


• C x+ ■ “ X+ Sy ' 


C =v fix , C . = v fix , 

y- y- * y+ y+ 


D = v6y/ 5x , D = v6y/ 5x , 

X™ XT 


Dy_ = y'Sx/'Sy , Dy + = vSx/5y . 


The term S( 4>) i j is the correction term and is implemented as 

S<*>i, 3 -1 [V»i-2,j - 2 +i-l,J + ♦i.J> 

+ V*l+2,j " 2 *i+l,j + *i,l ) 

+ C S<*l,J-2 - 2 *I,j-l + ♦i.J 5 

+ °N < j+2 ‘ 2 *l,j+l + ' 

Cy = max (0, C x _) , C E = max (0, -C^) , etc. 


where 


(18) 

(19) 

( 20 ) 
( 21 ) 


( 22 ) 


The value of S( <j>) is zero for the lines of cells i - 2, j - 2, i - IMAX-1 

and j = JMAX-1. P. . represents the finite-difference form of the pressure 

i > J 


term. 


Scheme 3 is implemented much like Scheme 1, but the values of the 
convective terms in A^, A^, etc. are doubled. Thus 

A w =■ max(0, 2C X _) + D x _ » 

A = max(0, -2C ) + D , etc. (23) 

E x+ x+ 


*For near-boundary cells the distances are half the cell dimensions for the 
velocity component parallel to the boundary. 
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The source term S( 4») . . is given by 

1 > J 

m 2 t<V - *l tj > 

+ °E ( *1+2,3 ' *1,3 > 
* C S <*1,3-2 - *1,3) 
+ C N <*1,3+2 - 


(24) 


Again, for the near-boundary cells, S( <J>) ± # j is set to zero and coefficients 
A^, Ag, etc., are reevaluated as in Equations (14) - (17). 


In contrast with the second order upwind differencing, the hybrid 
formulation uses a combination of first order upwind differencing and second 
order central differencing. In the hybrid formulation, central differencing 
is used on both convection and diffusion terms, but when the cell Reynolds 
number is greater than two, upwinding is used on the convection term and the 
diffusion is set to zero. The corresponding expressions for the 
coefficients are: 


A 


W 

A 


E 


A 


N 


max ( 

:|c 

l x- 

max ( 

:|c 

1 x+ 

max ( 

:|c 

1 y- 

max ( | C 

1 y+ 


, D ) + C 
x- x- 


• D x + } - C X + 


, D ) + C 

y- y- 


, D ) - C 
y+ y+ 


(25) 

(26) 

(27) 

(28) 


3. SOLUTION ALGORITHM 

In this study, the nonlinear finite-difference equations resulting from 
the second order upwind formulation are solved by a coupled multigrid 
algorithm. This algorithm was earlier used by the author in calculating the 
cavity flow, in conjunction with the hybrid differencing [25,26]. The 
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algorithm was observed to be rapidly convergent for large finite-difference 
grids (up to 514x514 grid nodes) and for very high Reynolds numbers (up to 
5000). 

The algorithm used here has two novel features. First, the momentum 
and continuity equations are solved simultaneously. Thus no pressure or 
pressure-correction equation is derived as is done in many earlier 
algorithms. Instead, the continuity equation is used in its primitive 
finite-difference form. Second, the coupled equations are solved 
iteratively by the multigrid technique [28]. The complete details of the 
present algorithm are given in references [25] and [26]. They are described 
here briefly. 

A. Solution Cycle 

In the present algorithm, the FAS-FMG (Full Approximation Storage - 
Full Multigrid) algorithm well suited for nonlinear problems is used. The 
FAS-FMG algorithm, flow charted in Figure 2, proceeds as follows. After a 
series of finite-difference grids is chosen, iterations are initiated on the 
coarsest grid (grid number 1). On this grid, the solution of the complete 
nonlinear problem is sought. The coefficients and the source terms in the 
finite-difference equations are evaluated with the existing values of 
velocities in store and the linear equations are then solved iteratively. 
The "relaxation procedure" used here is a Symmetrical Coupled Gauss-Siedel 
technique (SCGS). The SCGS was observed to provide good smoothing (i.e., 
convergence) rates for the hybrid scheme and is also used in this study. 
Basically, the SCGS scheme is a node by node solver; however, at each node, 
it updates all four velocities (i.e., on the four faces of the two- 
dimensional cell) and the pressure. The SCGS scheme solves five equations 
representing the complete transport characteristics at the node. In order 
to solve them, the equations are first written in terms of residuals and 
corrections at node (i,j). A matrix of four diagonal elements and two 
bordered line elements is then constructed from the five equations. Because 
of the bordered structure, the matrix can be easily inverted to give 
analytical expressions for the corrections. Due to the nonlinearity of the 
convection terms, the corrections are underrelaxed before adding to the flow 
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variables. The SCGS scheme is repeatedly applied at all nodes, updating the 
coefficients a n * Ag, etc., and the source terms within the iterative 
sweeps. Such sweeps are performed until the residuals have decreased below 
the prescribed tolerance level. The tolerance is applied to the summed 
averaged residual (suitably normalized) over the whole flow domain. 


After a converged solution on grid 1 is obtained, it is prolongated 
(extrapolated) to the next finest grid (grid 2). The solution on grid 2 is 
then sought. The coefficients are dynamically calculated inside the loop 
and at each node the five equations are solved. However, on this grid, the 
iterations are not carried until complete convergence. Instead, the rate of 
decrease of the averaged residual Is monitored and when this rate falls 
below a specified rate, i.e., when 



/ R 

P 


> n 


(29) 


(where subscript p here refers to iteration count), the iterations on this 
grid are temporarily halted. The residuals and the solution are then 
restricted to the next coarser grid. The coefficients are formulated on the 
coarser grid with the restricted values of the flow variables, and a 
solution to the coarse grid equations is sought. The equations solved on 
the coarser grid (say grid 1) are 


11 1 12 2 2 

L q 1 = F 1 + I* (r - Lq) , (30) 

1 2 2 2 

where restricts the residual (F - L q ) to grid 1. For the coarsest 
grid in the series of grids, Eq. (30) is solved exactly. 


The solution on grid 1 is now used to correct the solution on grid 2. 

At this stage, the solution on grid 1 is not prolongated to grid 2, rather 

the difference between the initial and final solution on grid 1 is 

2 

prolongated. Thus the correction to q (the superscripts here correspond to 
grid numbers and are not exponents) is 


2 2 2 
^new ^old + 1 




( 31 ) 
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The restrictions and prolongations between grids 1 and 2 are continued 

until an accurate solution on grid 2 is obtained. At this stage the 

solution on grid 2 is prolongated to grid 3 and the residuals are smoothed 

until the rate of smoothing falls below the threshold value (n). 

Restrictions are then made from grid 3 to grid 2 and the residuals are 

iterated on grid 2 until the smoothing is satisfactory. Otherwise the 

remaining residuals of grid 2 are transferred to grid 1 and annihilated 

-4 

completely (i.e. , to < 10 accuracy). As the solution sequence proceeds, 
finer and finer grids are considered and converged solutions are obtained on 
each grid. When a converged solution is obtained on the finest grid, M, the 
solution is terminated. 

_3 

The tolerance level on all grids except the coarsest is set to 10 
For the square cavity calculation, this represented a decrease of three 

orders of magnitude in the averaged residual. The tolerance on grid 1 is 

-4 

set to a smaller value (10 ) in order to obtain a nearly exact coarse grid 

solution (or correction). However, the accuracy level on intermediate grids 
during restriction from a finer grid is reset to a lower value. When an 
intermediate grid is reached by restricting the residuals from a finer grid, 

_3 

the solution on the intermediate grid is not obtained to 10 level of 
accuracy. Instead the tolerance is reset to 

£ h =Ae h+l’ (32) 

where A = 0. 2 and e^ + | represents the value of the residual for grid (h+1). 

B. Restriction and Prolongation Operators 

The restriction and prolongation procedures are somewhat dictated by 

the staggered mesh arrangement. Restriction is used for transferring fine 

grid values to a coarse grid, whereas prolongation is used for extrapolating 

a coarse grid correction to a fine grid. The two operators are denoted 
ti"“ 1 h 

by I, and I, , , respectively where h denotes one of the grids. A 
h h-1 

frequent restriction operator is injection, i.e., the coarse grid value is 
taken to be the local fine grid value. Thus 
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c f 

q i,j = q 2i-l,2j-l 


(33) 


where the superscripts c,f denote coarse and fine grid values. The 
injection operator is not applicable for the staggered mesh because of the 
different locations of variables on coarse and fine meshes (see Fig. 1). In 
the current study, the restrictions are made by averaging neighbor values. 
Let (ic,jc) and (if,jf) denote coarse and fine mesh indices, respectively. 
Also, let u i+1 / 2 j be referred to as u(i,j) and be referred to as 
u(i-l,j), etc. Then 


if = 2(ic) - 1 , and jf = 2(jc) - 1 , 

u C (ic, jc) = l/2[u f (if ,jf) + u f (if,jf-l)] , (34) 

v c (ic, j c) = l/2[v f (if ,jf) + v f (if-l,jf)] , (35) 

p c (ic,jc) = 1/4 [p f ( if ,jf) + p f (if,jf-l) , 

+ p f (if-l,jf) + p f (if-l,jf-l)] . (36) 


The prolongation relations are derived by a bilinear prolongation. For each 
coarse grid node, four fine grid values are calculated. The relations are 
given in Vanka [25). 


4. RESULTS 

In what follows, we shall present the results of calculations for the 
driven cavity flow situation using the combination of the second order 
upwind differencing and the coupled multigrid solution technique. A series 
of calculations employing several Reynolds numbers and finite-difference 
meshes has been completed. The calculations have been made with both 
Schemes 1 and 3 described above. The rates of convergence (i.e. decrease of 
residuals in the finite difference equations) are somewhat different for 
both schemes, but the converged solutions are the same. Therefore the 
converged solutions of only one of the schemes will be presented. The 
Reynolds numbers considered are 100, 400, and 600. For the last value, the 
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rate of convergence was quite slow and therefore larger Reynolds numbers 
were not considered. Several grids containing up to 80x80 finite difference 
cells were employed for each Reynolds number when convergence was 
possible. In order to test the relative stability of the second order 
upwind differencing vis-a-vis the hybrid formulation calculations were made 
also with the hybrid differencing. All calculations were initiated from 
zero velocity and pressure fields. 

A. Re = 100 


Calculations for Reynolds numbers of 100 were made with grids 
containing 10x10, 20x20, 40x40, and 80x80 internal cells. Further finer 
meshes were not considered because the 80x80 mesh provides grid independent 
results even with hybrid formulation. Figures 3-5 provide the rates of 
convergence of the three schemes for this Reynolds number. The iteration 
number corresponds to the fine grid and the residuals are the average summed 
normalized residuals in the momentum and continuity equation, defined as 

R = [I{(Rij) 2 + ( R I,j )2 + j ) 2 }/ (IMAX*JMAX*3) ] 1/2 . (37) 

The optimal relaxation factor for these calculations was 0.8 for all 
schemes. It is seen that the hybrid scheme converges fastest, followed by 
Scheme 1 and then Scheme 3. The CPU times for these calculations are given 
in Table 1. The CPU times, in general, are not proportional to the number 
of fine grid iterations because of the different number of iterations on the 
coarse grids. 


Table 1. CPU Times* (sec) for Re = 100 


Grid 

Scheme 1 

Scheme 3 

Hybrid 

10x10 

0.40 

0.59 

0.37 

20x20 

1.81 

2.26 

1.18 

40x40 

5.87 

7.89 

4.05 

80x80 

31.67 

30.01 

17.5 

*IBM 3033, 

F0RTHX, 0PT( 2) 

Compiler 
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The results of the converged fields are presented in Figures 6-8 for 
the different grids considered. For brevity, only the u-velocity profiles 
along the vertical centerline are presented. For comparison, the result of 
the 80x80 calculation is considered as the exact solution. From these 
figures, it is seen that for Re = 100, both the schemes give equally good 
accuracy with just 20x20 nodes. For the 10x10 grid, however, the second 
order upwind scheme is slightly more accurate. 

B. Re = 400 


The next Reynolds number considered is 400. Calculations for this 
Reynolds number were made with grids consisting of 10x10, 20x20, 40x40, and 
80x80 nodes. Another calculation with 160x160 nodes with the hybrid scheme 
is considered to be the exact solution. Figures 9-12 show the rates of 
convergence for the three schemes for the Reynolds number of 400. The 
hybrid scheme converged the fastest, in 14 iterations. The relaxation 
factor for the hybrid scheme was 0.8. For the second order upwind scheme, 
however, a more restrictive relaxation factor was necessary. On the coarser 
grids (10x10 and 20x20) a value of 0.6 was optimal. For finer meshes, it 
was necessary to reduce this value to 0.4. The number of iterations for the 
second order scheme was much larger than the number with the hybrid 
scheme. Thus, it is seen that the second order scheme is less stable at 
higher Reynolds number, a property that makes the scheme less attractive 
unless the additional accuracy gained more than compensates for the larger 
CPU time. 

Figures 13-16 show the centerline velocity with the hybrid and second 
order schemes for Re = 400, obtained with the various grids. The exact 
solution corresponds to the hybrid calculation with a 160x160 grid. This 
figure shows that both the second order upwind and the hybrid schemes give 
incorrect results even with 40x40 cells. However, with the hybrid scheme as 
the grid is refined further, the results get close to the exact solution. 

For the second order scheme, however, the solution does not approach 
the true solution in a monotonic way. For example, in some regions, the 
error in the 20x20 solution is more than that in the 10x10 solution. 
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Nevertheless, as the grid is refined to 80x80 nodes, the error decreases, 
but the solution still does not match the exact solution obtained with the 
160x160 grid and the hybrid differencing. This feature of second order 
upwind differencing is somewhat disturbing. The CPU times for the 
Re = 400 calculations are given in Table 2. 


Table 2. CPU Times* for Re = 400 


Grid 

Scheme 1 

Scheme 3 

Hybrid 

10x10 

0.66 

1.18 

0.50 

20x20 

5.375 

5.131 

1.93 

40x40 

25.01 

33.00 

7.10 

80x80 

90.50 

141.85 

22.81 


*IBM 3033, FORTHX, 0PT(2) Compiler 


C. Re = 600 


Calculations for Reynolds numbers of 600 were also made with grids 
consisting of 10x10, 20x20, 40x40, and 80x80. cells. For the hybrid scheme, 
the relaxation factor was 0.8, as before. However, the second order scheme 
was observed to require much more damping. It was necessary to lower the 
underrelaxation factor (equal values are used on both u and v velocities) to 
0.4 and 0.3 at finer grids. Figures 17-20 show the rates of convergence of 
the three schemes at this Reynolds number. Again, the hybrid is the fastest 
converging differencing scheme. Table 3 gives the required CPU times for 
these calculations. The results for Re = 600 are shown in Figures 21-24. 
Again the behavior is similar to that observed for Re = 400 but with 
slightly larger "overshoots" from the exact solution. 
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Table 3. CPU Times* for Re = 600 


Grid 

Scheme 1 

Scheme 3 

Hybrid 

10x10 

0.74 

1.56 

0.54 

20x20 

9.11 

8.43 

2.11 

40x40 

30.11 

30.19 

8.88 

80x80 

196.00 

192.52 

28.74 

*IBM 3033, 

FORTHX, 0PT(2) 

Compiler 



D. Higher Reynolds numbers 

Higher Reynolds numbers of 800 and 1000 were also considered, but it 
was observed that the above described trends, i.e., poor convergence and 
overshoots continue to be amplified with increase in Reynolds number. The 
hybrid scheme on the other hand was stable up to Re = 5000 with 160x160 
grids [25]. Further calculations with second order upwind scheme were 
therefore discontinued. 


5. DISCUSSION 

Our experiences with the second order upwind difference scheme combined 
with the coupled multigrid algorithm have been somewhat discouraging. At 
low Reynolds numbers we observe that the rate of convergence with second 
order upwind differencing is only slightly worse than with the hybrid scheme 
and the results are marginally better. At higher Reynolds numbers, the rate 
of convergence with the second order upwind scheme is significantly worse 
than the hybrid scheme. In addition, not much seems to have been gained in 
terms of improved accuracy. It is surprising that even with 80x80 nodes, 
the second order scheme has some error. Thus no clear benefits of using the 
higher order upwind scheme are observable. 

In contrast with our observations, Wilkes and Thompson [24] conclude 
that the second order upwind differencing is advantageous over the hybrid 
scheme (i.e., first order upwind scheme) even after discounting for the 
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additional time required for convergence. However, their studies were made 
in sudden expansion flows at low Reynolds numbers. Also in sudden expansion 
flows, the recirculation zone is smaller and of lower intensity than in the 
cavity situation. As a result, they observe improved performance of the 
second order upwind differencing over the first order (or hybrid) scheme. 

The conclusions by Shy and Correa [1] in a cavity flow are somewhat 
puzzling to us. They employed only 20x20 (equally spaced) cells but were 
able to calculate up to Reynolds number of 10,000. They employed the TEACH 
computer program, based on the SIMPLE technique. Our calculation procedure 
(which is observed to be very stable and rapidly convergent with the hybrid 
scheme) was not stable beyond Re = 800 unless prohibitively small relaxation 
factors are used. Also, the nature of the velocity profile presented by 
them (for Re = 1000) is different from our results at Re = 600. At this 
stage, the only difference we have identified between the two finite 
differencing procedures is their use of the hybrid formula at near-boundary 
cells versus our use of first order upwind formula. This appears to be a 
trivial difference and we do not feel this to be the reason behind the 
different observations. 

Although the present conclusions appear to discourage the use of second 
order upwind differencing, we believe that it may still be possible to use 
the scheme and take advantage of the second order accuracy if the precise 
cause of the "overshooting" is identified and remedied. We are of the 
opinion that some form of blending, with a first order scheme may be used in 
selected regions that have the potential for overshooting. A combination of 
first order and second order differencing may produce results superior to 
those of the individual schemes. We believe that further study in this 
direction will be fruitful. 
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Fig. 2. Solution Cycle 
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Fig. 5. Rates of Convergence, Re = 100, 40x40 Grid 
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Fig. 6. Comparison of Vertical Center-line u-velocity Profile for 

Re = 100, 10x10 Grid 
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Fig. 8. Comparison of Vertical Center-line u-velocity Profile for 

Re = 100, 40x40 Grid 





NORMALISED RESIDUAL NORMALISED RESIDUAL 


27 



Legend 

HYBRID 
SCHEME 1 
SCHEME 3 


Fig. 9. Rates of Convergence, Re = 400, 10x10 Grid 
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Fig. 10. Rates of Convergence, Re = 400, 20x20 Grid 
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Fig. 11. Rates of Convergence, Re = 400, 40x40 Grid 
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Fig. 12. Rates of Convergence, Re = 400, 80x80 Grid 
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Fig. 13. Comparison of Vertical Center-line u-velocity Profile for 

Re = 400, 10x10 Grid 



■j t i i i | 1 ! 1 1 1 1 f 1 i [ i i i i j i 1 1 r | 


0 0.2 0.4 0.6 0.8 1 

y 

Fig. 14. Comparison of Vertical Center-line u-velocity Profile for 

Re = 400, 20x20 Grid 




Fig. 16. Comparison of Vertical Center-line u-velocity Profile for 

Re = 400, 80x80 Grid 
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Fig. 19. Rates of Convergence, RE = 600, 40x40 Grid 
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Fig. 20. Rates of Convergence, Re = 600, 80x80 Grid 
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Fig. 22. Comparison of Vertical Center-Tine u-velocity Profile for 

Re = 600, 20x20 Grid 
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Fig. 23. Comparison of Vertical Center-line u-velocity Profile for 

Re = 600, 40x40 Grid 
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Fig. 24. Comparison of Vertical Center-line u-velocity Profile for 

Re = 600, 80x80 Grid 
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